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Abstract. Domain walls, optimal droplets and disorder chaos at zero temperature 
are studied numerically for the solid-on-soHd model on a random substrate. It is shown 
that the ensemble of random curves represented by the domain walls obeys Schramm's 
left passage formula with k = 4 whereas their fractal dimension is dg = 1.25, and 
therefore is noi described by "Stochastic-Loewner-Evolution" (SLE). Optimal droplets 
with a lateral size between L and 2L have the same fractal dimension as domain 
walls but an energy that saturates at a value of order 0(1) for L ^ oo such that 
arbitrarily large excitations exist which cost only a small amount of energy. Finally it 
is demonstrated that the sensitivity of the ground state to small changes of order 6 in 
the disorder is subtle: beyond a cross-over length scale Lg the correlations of the 

perturbed ground state with the unperturbed ground state, rescaled by the roughness, 
are suppressed and approach zero logarithmically. 
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1. Introduction 

Domain walls in disordered systems play an important role in understanding the stability 
of the ordered phase, the energetics of large scale excitations, the asymptotic dynamics 
in and out of equilibrium as well as the sensitivity to changes of external parameters. 
They have been studied quite intensively in the recent years for Ising spin glasses 
[Il[2l[3l[4l[5l[6l[7l[8],XY spin glasses [9], random field systems [lOl [IH [IS [l3] , random 
ferromagnets [lO], disordered elastic manifolds [ISl [IQ Ull [HI [HI ED] , and many others. 

Two domain wall properties are prominent: the first concerns energy and can be 
characterized by the scaling behavior of the domain wall energy with their lateral size, 
which gives rise to a first, sometimes universal exponent, the stiffness exponent 6. The 
second concerns geometry and gives rise to another, sometimes universal exponent, the 
fractal dimension ds, or, in case the domain is not fractal, a roughness exponent (. The 
interplay between energetics and geometry of the domain walls (i.e. between stiffness 
exponent and fractal dimension) determines how sensitive the system state is to changes 
of either external parameters like the temperature or a field, or internal parameters like 
small disorder variations. This sensitivity is often extreme in glassy systems and goes 
under the name of "chaos" [1]. 

Domain walls of glassy systems in two space dimensions represent fractal curves 
in the plane and the question arises, whether they fall into the general classification 
scheme for ensembles of random curves described by Stochastic Loewner Evolution 
(SLE) [2H [22]. Recently indications were found that domain walls in 2d spin glasses 
(at zero temperature) are indeed described by SLE [23], [24], at least for a Gaussian 
distribution of the bonds, but apparently not for binary couplings [25]. Also the domain 
walls in the random-bond Potts model at the critical point (i.e at finite temperature) 
were found to be numerically consistent with SLE [26]. It appears natural to ask, 
whether the domain walls in other two-dimensional disordered systems are potential 
candidates for a description by SLE. 

In this paper we study domain walls and chaos at zero temperature in the solid- 
on-solid (SOS) model on a disordered substrate. This is a numerically convenient 
representation of a two-dimensional elastic medium, with scalar displacement field, 
interacting with quenched periodic disorder. It has been studied to describe various 
physical situations ranging from vortex lattices in superconductors to incommensurate 
charge density waves and crystal growth on a disordered substrate [27l [28l [29] . Here we 
focus on three questions: 1) are domain walls in this model described by SLE, 2) what 
is the relation between size and energy of optimal excitations (droplets) in this model, 
3) does disorder chaos exist in the ground state of this model ? After a brief summary 
of what is already known about the model and a description of the numerical method by 
which we compute the ground state and the domain walls these three issues are studied 
in separate sections. The paper ends with a discussion of the results obtained. 
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1.1. Model 

We consider the solid-on-solid model on a disordered substrate defined by the 
Hamiltonian 

H = ^(/ij - hjY , hi = ni + di , (1) 

with i = {xi,yi) G Z^. In Eq. the height variables rii {i = 1,...,N) take on 
integer values rii = 0, ±1, ±2, . . . and the offsets di are independent quenched random 
variables uniformly distributed between and 1. The sum is over all nearest neighbor 
pairs (ij) of a rectangular lattice of size x Ly (L^ = Ly = L if not stated otherwise). 
The boundary conditions will be specified below in the context of domain walls. The 
Hamiltonian in Eq. (P) describes a discrete model of a two-dimensional elastic medium 
in a disordered environment. In the continuum limit, it is described by a sine-Gordon 
model with random phase shifts (and in the absence of vortices), the so called Cardy- 
Ostlund model [30] . 

Hco = J dh{Vuir)f - A cos(27r[M(r) - d{r)]) , (2) 

with a continuous scalar displacement field u{r) G (— oo, +oo) and quenched random 
variables (i(r) G [0,1]. Discretizing the integral and performing the infinite strong 
coupling limit A — > oo one recovers 

It is well known that this model ([H [2]) displays a transition between a high 
temperature phase, T > Tg = 2/7r where the disorder is irrelevant and a low temperature 
phase below Tg, dominated by the disorder. The high-temperature phase T > Tg is 
characterized by a logarithmic thermal roughness 

C{r) = {{h,-h,+ry)-^T\ogr, (3) 

where (. . .) denotes the thermal average and the average over the quenched disorder. 
The low-temperature or glassy phase is instead "superrough", characterized by an 
asymptotically stronger (log-square) increase of C(r): 

C(r) ~ c(T) ■ log2 r + 0(log(r)) , (4) 

which means C = 0, as expected for a random periodic system. Close to T^, a Coulomb 
Gas renormalization group (RG) analysis to lowest order gives c(T) ^ (1 — T/Tg)"^ /2tt'^ 
[29], in rather good agreement with numerical simulations |31| . At T = 0, numerical 
simulations give the estimate c(T = 0) ^ 0.5/(27?)^ ^ 0.012 ^ |32]. While earlier 
studies, based on "nearly conformal" field theory |33|, claimed an exact result for A(T), 
predicting A{T = 0) = 0, in clear contradiction with numerics, a more recent approach 
based on PRC, incorporating non analytic operators predicts a non-zero A(T = 0) which 
compares reasonably with numerics |34| . 

For free or periodic boundary conditions, the Hamiltonians ([I]) and ([2]) have a 
discrete symmetry, the energy is invariant under a global height (displacement) shift 
Hi ^ Hi + An {u{r) u{r) + An), where An is an arbitrary integer. This symmetry 
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Figure 1. Left: Ground state n'' of a 200 x 200 system with the boundary sites 
(indicate in green) fixed to rii = 0. The different height values rii are grey-coded 
(dark = low values, bright = high values). Middle: Ground state configuration of the 
same system as to the left with the upper half of the boundary sites (indicated in red) 
fixed to ni = 1 and the lower half (indicated in green) to = 0. Right: Difference 
plot between the Left and Middle plots: in the lower white region the ground state 
configuration is identical to the corresponding sites in the left figure, whereas in the 
upper grey region they differ by exactly An = 1 from the corresponding site in the 
middle panel. The border between the white and the grey region is a domain wall, 
representing a step in the height profile of the ground state. 



will not be broken in the low temperature phase of the infinite system and true long- 
range order at T < T^, is absent, i.e. {hi) = 0. Concomitantly the model (P), with free 
or periodic boundary conditions, has infinitely many ground states, which differ by a 
global shift An G {±1, ±2, . . .}. 

1.2. Domain walls 

By an appropriate choice of boundary conditions one can force a domain wall into the 
system, which is most easily visualized at T = (c.f. Fig. [J]): consider the square 
geometry and fix the values of the boundary variables to rij = 0. This yields a unique 
ground state configuration If one fixes the boundary variables to nj = +1, the 
corresponding ground state would be n'° = n° + 1. A domain wall inducing boundary 
condition is one, in which the lower half of the boundary values are fixed to = and 
the upper half to n° = 1. The ground state of this set-up is then fi^ = n° in some, 
mainly the lower region of the system, and = in the rest - both region separated 
by a domain wall of non-trivial shape. 

It turns out that these domain walls are fractal [HI [H], which means that their 
lengths Zdw scales with linear system size as 

Uh ~ , (5) 

with dg > d — 1 = 1. The numerical estimate for dg is dg = 1.27 ± 0.02 |jl8j. Such 
a fractal scaling of zero-T domain walls is also found for spin glasses, in the 2d EA 
model with Gaussian couplings it is dg^sG = 1-27 ± 0.01, and with binary couplings it 
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is ds^sGB = 1-33 ± 0.01. On the other hand zero-T domain walls in disordered Ising or 
Potts ferromagnets are rough (i.e. are characterized by algebraic correlations) but not 
fractal. 

The energy for such a domain wall, given by the difference between the energy of 
the ground state of the system with the domain wall inducing boundary conditions and 
the one with homogeneous boundary conditions, increases with L logarithmically 

A^~logL. (6) 

This result, which was obtained by numerical simulations [16], is consistent with the 
usual scaling relation AE ~ together with the exact result 6 = d — 2 + 2( = 
(thanks to statistical tilt symmetry [35]). This logarithmic behavior is characteristic of 
a marginal glass phase, described by a line of fixed point indexed by temperature (which 
is here marginal in the RG sense). For comparison the stiffness exponent in 2d (3d) spin 
glasses is OsG2d = — 0.28±0.01 and OsG2d = 0.3±0.1 (and thus characterized by a T = 
fixed point), whereas for disordered Ising or Potts ferromagnets AE ~ L'^^^. 



1.3. Method 

The ground states of i.e. the configuration n° = (n5,...,n^) with the lowest 
value for the energy /7[n°] for a given disorder configuration d = (rfi, . . . , ^at), can 
be computed very efficiently using a minimum-cost-flow-algorithm [36l UHl [37] . For the 
speciflc details in which domain walls are induced in the ground state it is useful to 
recapitulate the mapping onto a minimum-cost-flow problem. 

After introducing the height-differences n*^ = rii — uj (integer) and d*j = dj — di 
(e [— 1,+1]) along the links k = {i,j) on the dual lattice G* one obtains a cost (or 
energy) function that lives on the dual lattice 

k 

The conflgurations = (nj, . . . ,n%j), where M is the number of links (or bonds) of 
the original lattice, constitute a "flow" on the graph G*. Suppose the original model 
([I]) has free boundary conditions. Then the sum of the height differences along any 
directed cycle in the original lattice vanishes. Therefore the divergence of n* vanishes 
at all sites i: 

(V • n*), = , (8) 

which means that the flow n* on G*, in order to give rise to a height field n on the 
original lattice G has to be divergence-less, i.e. without sources or sinks. The problem 
of determining the ground state n of ([I]) is thus equivalent to find the flow n* with the 
minimum cost 1^ under the mass-balance constraint ((HI) - i.e. a minimum cost flow 
problem, for which there exist very powerful algorithms [36l [TBI [37]. 

Enforcing one domain wall, or step of height one, into the ground state of ^ 
by appropriate boundary condition is then equivalent to modify the constraint ([8]) at 
exactly two sites, the start and end point of the domain wall (see Fig. [2h-d). As an 
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Figure 2. Different geometries and constraints on domain walls considered here: a) 
boundary conditions inducing a step / domain wall as in Fig.[TJ b) Boundary conditions 
inducing a step / domain wall running diagonally from one corner of the rectangular 
lattice to the opposite one. c) Boundary conditions for a circular domain inducing a 
boundary along the equator with two different orientations of the underlying lattice, 
d) Boundary conditions for a half circle domain and a domain wall with one fixed end 
at the origin and a free end on the outer half circle. 

example consider the case in which one wants the domain wall to start at the point 
{x,y) = (l,L/2) and end at {x,y) = (L, L/2) of a square lattice. Then one chooses the 
boundary conditions for Ui as follows (c.f. Fig. [2] a): one fixes the values for rii at the 
lower half of the boundary (i.e. at z = (x, 1) for a; = 1, . . . , L and i = (1, y) and (L, y) 
ioi y = 1, . . . , L/2) to nj = 0, and the values for rii at the upper half of the boundary 
(i.e. at z = (x, L) for X = 1, . . . , L and at i = (1, y) and (L, y) for y = L/2 + 1, . . . , L) to 
rzj = 1. Translating these boundary conditions for the height variables n into constraints 
for the flow variables n* one immediately sees that at the point {x,y) = (l,L/2) and 
(x, y) = (L, L/2), where the step in the height profile starts and terminates, respectively, 
the constraint (l8|) is modified into 

(V-n*)(i,i/2) = +l, {V-n*)^L,L/2) = -l. (9) 

In other words: the induced step sends a unit of fiow from the starting point of the 
domain wall, which is the a source of unit strength, across the sample to the end point, 
the sink, along an optimal (minimum cost/energy) path. In what follows we identify 
domain walls immediately with the optimal path for the extra fiow unit defined by the 
modified mass balance constraints 

With the help of this concept one can then also consider situations in which the 
starting point of the domain wall is fixed but the ending is only forced to be on a specific 
region of the boundary, opposing the starting point (see Fig. [2li). Suppose one wants 
he domain wall to start at is = (l,L/2), and terminate somewhere on the opposing 
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boundary if = {L,y) with ?/ G {1, . . . ,L}. Then one introduces an extra node into the 
dual graph G*, denoted as the target node, connects it with bonds of zero cost to all 
sites on the terminal boundary, and assigns a sink strength —1 to it. The source node is 
the one closest to is in the dual graph and has source strength +1. The minimum cost 
flow of this arrangement is then the desired ground state configuration with a domain 
wall starting at ig and ending somewhere on the opposite boundary. 



2. Schramm-Loewner evolution (SLE) 

Since the domain walls as defined above represent fractal curves embedded in a two- 
dimensional space the question arises whether they fall into the classification scheme 
of Schramm-Loewner evolution (SLE) like loop-erased random walks, percolation hulls, 
and domain walls at phase transitions in 2d in the scaling limit [22l [2T]. The necessary 
(and sufficient) condition for a set of random curves connecting two points on the 
boundary D of a domain to be described by SLE are 1) the measure for these random 
curves has to fulfill a Markov property, 2) the measure has to be invariant under 
conformal mappings of the domain. Recently it was suggested that also domain walls in 
2d spin glasses can be described by SLE [23l [24] although neither conformal invariance 
is fulfilled for each individual disorder realization nor the Markov property after disorder 
averaging. However, conformal invariance might hold for the disorder averaged model 
and the Markov property might be fulfilled almost always in a statistical sense. 

A single parameter k parameterizes all SLEs and it is related to the fractal 
dimension of the curve via 

4 = 1 + k/8 . (10) 

The parameter k is the diffusion coefficient of the Brownian motion that underlies the 
SLE and generates via a random sequence of simple conformal maps the fractal curves. 
In addition to the fractal dimension it determines various other geometric and statistical 
properties of the SLE curves. One of them is for instance the probability that a curve in 
the upper half plane EI generated by SLE will pass to the left of a given point z = x + iy 
is given by Schramm's "left passage formula" [38] 



1 Tii) „ (i i:i 



2^ 



oo 



where 2-^1 is the hyper-geometric function 2-F1 (01,02; 6; 2) = Yl ^rfaiT^ ^nlT^ ^ • 

k=0 

Since the probability depends just on the ratio between Re{z) and Im(z), it is sometimes 
useful to replace this ratio by a function of an angle. We decided to use: tan (0) = x/y. 
So </) g] — |; |[ is the angle at the origin between the imaginary axis and z. This formula 
holds for the domain of the SLE being the upper half plane with the start point of the 
curve being identical to the origin of the coordinate system and the end point at infinity. 

A standard check whether an ensemble of random curves is a potential candidate 
for SLE therefore is to simultaneously determine their fractal dimension dg and the 
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Figure 3. Average length of domain walls spanning the system from one end to the 
opposite one as a function of the system size L in a log- log plot. The straight lines are 
least square fits to ^ and yield the estimate for the fractal dimension ds — 1.25±0.01. 



left passage probabilities, and to test whether the latter fit (fTTll with k = 8{ds — 1) 
[231 El EH]. 

In Fig. [3] we show our data for the average length of a domain wall in the 
configuration depicted in Fig. [2] b, i.e. starting at the point ig = (1, 1) and ending 
at it = (L, L) in a L X L geometry and at it = {L, 2L) in a L x 2L geometry. Least 
square fits to the scaling law (I5|) yield the estimate dg = 1.25 ± 0.01 which agrees with 
the value found in [T8| . 

This value for the fractal dimension would imply k = 2.00 ± 0.08 if the domain 
walls are described by SLE. Next we determined for different points (x, y) of the lattice 
the frequency that a domain wall passes to the left of it, yielding a probability p{x, y), 
which we compared with Schramm's left passage formula P^{x,y) for fixed k. For the 
circle domain and the choice of the start and end points of the domain wall as shown in 
Fig. [2] c, the formula (fTB is modified: 

Let E be the unit circle in the complex plane. The Cayley function g : K ^ 
H, z I— > i^rf maps the unit circle conformally into the upper half plane H. Furthermore 
g{—l) = 0, g{l) = oo, thus curves in E starting at z = — 1 and ending at 2; = 1 are 
conformally mapped on curves in H. For the latter, if the curves are described by SLE, 
Schramm's formula ( fTTll holds, so that for the former, the modified formula 

P (z)-\^^&— F 4 3 /Regiz)y\ Regjz) 
""'^ ^ " 2 + 0Fr (^) '"^^ [r 2' [lmg{z)) ) lmg{zy ^''^ 

holds. Hence we filled the unit circle with finer and finer grids G approximating better 
and better the continuum limit for the curves in E that we want to check for SLE. For 
this we define the function 



fG{^^)= J2 [p{x,y)-P>^A^,y)]' (13) 

that measures the cumulative squared deviation of the probabilities p{x, y) from the 
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Figure 4. Left: The cumulative squared deviation fcii^) of the computed left 
passage probabilities p{x,y) from the values PK,.g{x,y) given by (fT2l) as a function 
of K. The underlying lattice geometry is the circle as sketched in Fig. [2fc., the 
corresponding conformal map g{z) entering (fT2|) is given in the text. The minimum is 
at K = 4.00 ± 0.01 with a squared difference per grid point of about 2 • 10^^. Right: 
Absolute difference between the calculated left passage probability p{x, y) and the SLE 
expectation P^^gix^y) (fT2| for k = 4 as a function of the 2d lattice coordinates {x,y). 
For the whole unit circle the deviation is almost everywhere smaller than 1%. Note 
that the domain wall is fixed at (—1,0) and (1,0), where the largest deviations occur. 




3 3.5 4 4.5 5 



Figure 5. Cumulative deviation fciK) as function of k and absolute difference between 
p{x,y) and Pi^,g(x,y) for k = 4 as in Fig. |4] but for the square geometry as depicted 
in Fig. [2}3. Note that the domain wall is fixed at (0,0) and {L,L), where the largest 
deviations occur. 

SLE- value for given k. The result is shown in Fig. [H The minimal deviation of the 
data from the expected SLE result is at k = 4.00 ± 0.01 which clearly differs from the 
value K = 2.00 ± 0.08 that one would expect from the fractal dimension if the domain 
walls would be described by SLE. 

Next we varied the geometry and the domain wall constraints and studied the 
case depicted in Fig. [2]b, i.e. a quadratic domain D with corners at 0, p + ip, 
ip {p real and positive). The function g : D ^ M. : z —p{z]S) with 5* = 
{2nip + i ■ 2n2p\ rii, n2 G defines a conformal map from D into H, where p{z] S) = 
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Figure 6. Deviation of the left passage probability p((/)) from the SLE expectation 
Pk{4>) ifTTji for K = 4 in the half circle geometry depicted in Fig. [2]d (p = 510, q = 
255, R — 500), i.e. one free end of the domain wall. 



is the Weierstral?. p-function. Furthermore g{p+ip) = 0, g{0) 

w£S\{o} 

oo, thus curves in D starting at z = and ending at z = p + ip are conformally mapped 
on curves in H starting at the origin and extending to infinity. Using g in (fT2l) we 
determined k in the same way it was done in the case of the unit circle. The result is 
shown in Fig. [5] and yields n = 4.00 ± 0.01. 

We also compared p{x, y) directly with Pf^=i{x, y) for curves in the upper half plane 
EI starting at the origin. For this we considered the geometry depicted in Fig. [2li, in 
which the curves starting at the origin can end everywhere on the half circle. We checked 
the prediction (fTTl) for different radii R and angles $ = arctan G ]— f ; f [• The 
result is shown in Fig. [6l where the deviations from Pk{4>) ([TTll for k = 4 are everywhere 
less than 2% for the size L = 500 and radii shown in Fig. [6l We also observe that the 
deviation systematically decreases with the system size L for fixed ratio R/L indicating 
vanishing deviations -P($) — -Pk=4($) in the limit L ^ oo. 

For all the geometries that we studied we found that the data can be nicely fitted 
with K = 4, but the fractal dimensions of the domain walls is unchanged by the different 
set-ups: ds = 1.25 ± 0.01. The conclusion is that domain walls in the SOS model on a 
disordered substrate are not described by SLE. The question arises which condition for 
SLE is actually violated. The fact that for all geometries the left passage probability 
we studied is well described by a common expression, Schramm's formula containing 
the conformal map of the half plane to the specific geometry under consideration, would 
actually hint at conformal invariance. But obviously this does not exclude the possibility 
that a breaking of conformal invariance manifests itself in other quantities, or even other 
geometries. We did not attempt to check the domain Markov property, as was tried 
in El. 
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3. Excitations 

In this section we study large scale excitations that cost a minimum amount of energy, 
also denoted as droplets [2], and the size dependence of their energy. According to the 
usual arguments in droplet scaling theory [2] this excitation energy is expected to scale 
in the same way as a domain wall of lateral size L, i.e. like AE ~ with 9 the stiffness 
exponent, which for the system we consider is = 0. The domain wall energy scales as 
Ai^DW ~ logi^ c.f. ([6]), and if it is correct that all excitations of scale L display the same 
behavior the question arises: how can thermal fuctuations destabilize the ground state 
such that the present model is indeed characterized by a line of fixed points, indexed by 
temperature (and thus different from a T = fixed point like in spin glasses with d > 21 
A logL scaling of excitations still implies that larger and larger excitations, necessary 
to occur to decorrelate the system's state from the ground state, need more and more 
energy and are therefore more and more unlikely, i.e. occurring with a probability that 
decays as exp(-AE/T) ~ L'"/^ at T > 0. 

The solution to this problem lies in the freedom that domain walls of excitations 
of scale L have to optimize their energy (c.f. for a similar discussion in the context of 
optimized dislocations in [IS]). Finding the optimal excitation of a given scale requires 
a higher effort than searching a domain wall with given start and endpoints (as we have 
seen above, there are even efficient ways to optimize the positions of these endpoints). 
There is a simple way to induce excitations of an arbitrary scale, as first proposed in the 
context of spin glasses [4] and used again in [8]: we could compute the ground state for 
fixed boundary conditions as depicted in the left panel of Fig.[T], this yields n°. Then we 
choose a central site (say i = {x, y) = {L/2, fix it to rij = + 1, fix the boundary 

as before and compute the ground state again, giving n°. This will differ from n° only 
by a compact cluster C that contains the site i and which has n° = n° + 1. This is then 
in fact an optimal or droplet excitation, but its size V = #{z G C} can vary from 1 to 
L^. Droplets of a fixed size cannot be generated in this way. 

The SOS model on a disordered substrate actually allows for an efficient search for 
optimal excitations of a given scale, as we describe in the following: 

A droplet excitations in the SOS model is a simply connected cluster C of sites 
whose height values are all increased (or all decreased) by one as compared to the ground 
state: = ?t,° + 1 Vz G C, see Fig. [7l Pictorially it is, in the height representation, an 
extra-mountain (valley) of height +1 (—1). Its boundary is a directed cycle in the dual 
lattice, and remembering the mapping to the minimum-cost-fiow problem (JTHH]), adding 
a cycle to a feasible solution maintains the divergence-free constraint ([8]). Hence, the 
transition from one state to the other is determined by such a cycle. Thus first one has 
to compute the ground state of a given disorder realization. Then in order to find a 
droplet excitation of this ground state that has a given lateral size one can for instance 
force this extra-cycle to run within an annulus of inner radius L/4 and outer radius 
3L/4 (i.e it's average diameter is L/2). This can be achieved by simply removing all 
sites / bonds outside this annulus and then computing the optimal cycle within this 




Figure 7. a) Sketch of a droplet excitation representing a simply connected cluster C 
of sites i whose height values are increased by one as compared to the ground state: 
rii = n-* + 1 Vi e C. b) Directed cycle in the dual lattice corresponding to the cluster 
in a) representing the step in height profile induced by raising the cluster C by one 
height unit, c) Sketch of an s — t cut of the nodes of the original lattice graph such 
that those sites connected to the external node s are forced to be in the set S and 
those connected to t in the set T (see text). 



modified graph, the cost for which depends on the ground state configuration and the 
substrate heights. One assigns costs to each directed edge that corresponds to the 
energy cost for increasing the height difference between its left and right side by one 
unit [36l [HI [37]. Note that in practice one may use the reduced costs emerging from 
the ground state. That is, if and only if a feasible fiow has minimum energy, then there 
are non-negative reduced costs such that the costs of each cycle remains unchanged. If 
the successive shortest path algorithm is used to solve the minimum-cost fiow-problem 
when computing the ground state n°, then no extra work is necessary to compute the 
desired non-negative reduced costs. Hence, Dijkstra's shortest path algorithm can be 
used to find the shortest directed cycle around the annulus, i.e. the one separating the 
inner and the outer ring of the annulus. To this end, the annulus is cut from the outer 
to the inner ring, i.e. the corresponding edges are removed from the graph, to prevent 
the shortest path from short-cutting. For each of these removed edges, the shortest path 
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Figure 8. Examples of optimal excitations of scale L, whose boundary (domain wall) 
is forced to lie in the interior of the area indicated in green (an annulus in square 
geometry). From the top row to the bottom it is L = 32, 64, 128 and 256. Pictures 
are scaled to have the same size. 



from its head to its tail is computed, where only the the remaining edges are used. The 
obtained shortest paths are completed with the corresponding directed edges to form 
cycles around the annulus. This procedure has to be repeated for all (or a representative 
number of) positions of the annulus within the original lattice (in practice one fixes the 
annulus and shifts the disorder configurations, wrapping it around a toroidal geometry). 

The procedure of finding the optimal cycle in a given annulus can be simplified 
by observing that the droplet excitation in the height representation corresponds to an 
s-t-cut of the underlying graph [37] in such a way that one forces all nodes of the 
inner circle of the annulus to belong to S and all nodes belonging to the outer ring of 
the annulus to T, see Fig. [7| c. The minimum s-t-cut with respect to the non-negative 
reduced edge costs of the ground state n° is then exactly the boundary of the optimal 
excitation (or the optimal cycle) one is searching for a given annulus arrangement. 

I An s-t-cut is a partition of the nodes of a graph Q into two disjoint sets S and T — Q /T, such that 
s e S" and t e T. 
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Figure 9. Disorder averaged energy of the optimal excitations of scale L as a function 
of L - left in linear scale, right in a log-log scale. 



According to the famous Min-Cut-Max-Flow theorem one can compute the minimum 
s-t cut in polynomial time by solving the associated maximum-flow problem [TBI [37] . 
We have implemented this procedure and show for illustration a number of examples in 
Fig. El 

Fig. [9] shows our result for the disorder averaged energy of the optimal excitations 
of scale L that we obtain with the procedure described above. This represents an upper 
bound for the optimal excitations of scale L since the annulus arrangement does not 
include all possible excitations of scale L. As one can see this bound saturates at a flnite 
energy of order 0{1) in the limit L — > oo. Consequently arbitrarily large excitations 
exist that cost only a small amount of energy, which renders the ground state unstable 
at an non-vanishing temperature. 

Fig. [ini shows the distribution of optimal excitation energies for different values of 
L. As can be seen the distribution is nearly Gaussian with a flnite width, i.e. it does 
not display long, e.g. algebraically decaying, tails, which implies that the average is 
representative for almost all disorder conflgurations. It can also be seen that the whole 
probability distribution Pl{AE) becomes independent of the length scale L for large 
L. This is an important observation since droplet scaling theory ^ one would expect 
Pl{AE) ~ l~^p{AE/L^). For 6* = 0, as it is the case here, it is not a priori clear whether 
this implies Pl{AE) ~ (In L)~^p{AE / In L) , i.e. a scaling with the average domain wall 
energy InL, or Pl{AE) ~ p{AE), i.e. droplet size independence. Fig. [TO] shows that the 
latter is correct. 

This has important implications for the scaling of the average droplet energy 
within a system of lateral size L as determined in [201. There the average energy of 
droplets of size / < L, (Ai5)™", i.e. without a lower bound, was estimated to behave 
like (A-E)™™ ~ InL. In order to make contact with our result one should note that 
this energy is the minimum among droplets of the kind we determined here, with size 
/ G [L/2,L], I e [L/4, L/2], / e [L/8,L/4], . . ., i.e a minmum of approximately InL 
random numbers. Only if the probability distribution of these energies is a) independent 
and b) identically distributed, their minimum goes like the inverse of their number, i.e. 
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{AE)f'' = min{(AE)i/2i, {AE)l/2^, . . . , {AE)L/2k} ~ l/k with ~ InL. According to 
our result depicted in Fig. [10] the assumption b, which is impUcit in the reasoning of 
[20] , is indeed fulfilled. 

Finally we note that also determined the fractal dimension of the boundary of 
the optimal excitations and found that it is identical with the fractal dimension of the 
domain walls with fixed start and endpoints: dg = 1.25 ± 0.01. 

4. Disorder Chaos 

In this section, we study the sensibility of the ground state to a small change of the 
quenched disorder configuration. To this purpose we generate a configuration of random 
offsets d} and compute the associated ground state hj. Then we slightly perturb this 
configuration of random offsets df = di + Sei with 5 <C 1 and where e^'s are independent 
and identically distributed Gaussian variables of unit variance and we compute the 
associated ground state hj. The question we ask is : how different are these two 
configurations of the systems h] and /if ? 

Such questions first arose in the context of spin glasses [T], where it was proposed 
that disorder induced glass phases may exhibit "static chaos", i.e. extreme sensitivity 
to such small modifications of external parameters (like disorder considered here or 
temperature). Such small perturbations are argued to decorrelate the system beyond 
the so called overlap length Ls which diverges for small 5 as ~ with a the 

chaos exponent. As an example let us first consider the case of an Ising spin-glass 
with a continuous Gaussian distribution of random exchange interactions, of width 
J. The system has two ground states related by a global spin reversal. Within the 
phenomenological droplet theory [D, [HJ [H] , a low-lying energy excitations of the system 
involves an overturned droplet of linear size L and costs an energy JL^ . If we now add 
a small random bond perturbation, say Gaussian of width 5J , the excess energy of 
a droplet is modified. For such spin-glass system, this energy comes only from the 
bonds which are at the surface of the droplet. Their contribution is thus the sum of 
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L'^" independent random variables of width 6 J, where dg is the fractal dimension of 
the droplet : it is thus of order ±51''"'^'^. Therefore the ground state is unstable to 
the perturbation on length scales L such that dJL''"'/'^ > JL^, i.e. L > Ls where 
Ls ~ (5~^/"sG -v^^ith asG = ds/2 — 9. One thus sees that, for spin glasses, disorder chaos 
is closely related to the (geometrical) properties of the domain walls. 

The situation is rather different for elastic systems in a random potential as 
considered here. Here we consider the SOS model on a disordered substrate defined 
in Eq. ([I]) as 

H = - hjf , hi = ni + di, (14) 

with i = {xi,yi) G Z^. In Eq. (fT4l ) the height variables Ui (i = 1,...,N) take on 
integer values Ui = 0, ±1, ±2, . . . and the offsets di are independent quenched random 
variables uniformly distributed between and 1. This system (with free of periodic 
boundary conditions) has infinitely many ground states which differ by a global shift 
An G {±1, ±2, ...}. Again, within the droplet argument [HI [42l [43] a low-lying energy 
excitation of the system involves a droplet of size L where the height field is shifted by 
unity, say An = 1 and it costs an energy . But now if we had a small perturbation 
(if = d} + 6ei, the excess energy of this droplet comes from the bulk of the droplet, which 
experiences a random force field. This contribution is thus the sum of L'^ random 
variables of width 6 and therefore in this case the ground state is unstable to the 
perturbation on length scales L > Ls where Ls ~ 5^^^" with a = d/2 — 9. Here 
d = 2 and 9 = and thus one expects Ls oc 6~^. At variance with spin-glasses discussed 
above, one thus sees that for disordered elastic systems, disorder chaos is not directly 
related to the properties of domain walls. 

For the present model (fT4ll . disorder chaos was demonstrated analytically at finite 
T near the glass transition Tg using a Coulomb Gas Renormalization Group [44]. At 
T = 0, some indications of disorder chaos were also found numerically in Ref. [16] where 
global correlations between hj and /if where studied through x((5) = J^ii^l " ■ 
this paper, following Ref. [43], we characterize the local correlations between these two 
configurations by the correlation function Cij{v) with r = (x, y) (with i,j = 1,2) 

a,(r) = {K-hl_,^{hi-hi_,^ , (15) 

where k + r = {xk + x,yk + y) and for a rotationnaly invariant system considered here 
one has Cjj(r) = Cij{r) with r = |r|. In the following we will talk about intralayer 
correlations for Cair) and infer/tt?/er correlations for Ci-/,j{r). Equivalently, one can also 
study such correlations (fTSl) in Fourier space and define Sij{c\) with q = {qx,<iy) as: 

k 

where q ■ = qxXk + qyyk- For a rotationnaly invariant system considered here, one has 
Sij{c\} = Sij{q) where q= |q|. 

Disorder chaos at T = in generic disordered elastic systems in dimension d was 
recently studied analytically using the Functional Renormalization Group (FRG) [43] . 
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At one loop order in a dimensional expansion in c? = 4 — e it was found that for short 
range disorder (like random bond problems) and random periodic systems in dimension 
d > 2 (including one component Bragg-Glass), one has [43] 

f c^* X -C 1 

Cu(r) = r2^$(5r") with $(x) ~ <^ _ ' ' (17) 

[ X ^ , X > 1 , 

with c^^ a constant, ( the roughness exponent and where n is the decorrelation exponent. 
Translated into Fourier space, this yields 

S^2{q) = Lf-^'^^^iqLs) with ^^(x) \ ^[^^ ' f ^ ' (18) 

[ X ^ , X > 1 , 

with Ls ~ 5^^/" and where the behavior for large x is then such that the dependence 
on Ls cancels in this limit, as it should. 

The two-dimensional disordered SOS model we are considering here (fT4l) 
corresponds precisely to the marginal case d = 2 (with C = and thus 6 = 0) where, 
as discussed in Ref. [43l |45], the analysis yielding the result in Eq. ffTTI) ceases to be 
valid. Indeed in that case non local terms, irrelevant in d > 2 are generated under 
coarse-graining and these additional terms have to be handled with care at T = 0. One 
thus considers the Hamiltonian associated to the two copies of the system parametrized 
by the scalar fields = M*(r): 

^(Vr^T + V,{U\ r) - fl\T) ■ VrU' , (19) 

1=1,2"' L 

where /x* = {jJi^xil^^) ^'^^ two-component random tilt fields, which are generated upon 
renormalization. While they are irrelevant in d > 2 they become relevant in d = 2 where 
they play a crucial role. These two copies v} , v? are thus independent (fTQ]) but they feel 
two mutually correlated random potentials 

Viiu, T)Vj{u', Y') = Rijiu - u')5'^{t - y') , (20) 

= ^'/pp"^'(r - r') , (21) 
where z = 1, 2 is the index of the copy and p = x, ?/ is a spatial index. This leads to the 
replicated Hamiltonian, Z'^ = exp (— ifrep): 



-f^2copies — ^ 2^ I d Y 



H.e, = ^ fd'rJ2t. - ^ E E / d'mM - ^^22) 

1=1 a=l i,j=l a,b=l 

where we used the notation = u\ — u\. In the "bare" model, one has Gij{u) = aij. 

Close to the transition T < Tg, this model (l22l) was studied using Wilson RG 
analysis by varying the short scale momentum cutoff = Ae~^, with i the log-scale. 
It was shown, using Coulomb Gas technique at lowest order, that Gu oc a£. This 
leads to the correlation function in Fourier space Sii{q) oc cri/q"^, which yields, setting 
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£ = log(l/g), to the log^ (r) behavior of the intralayer correlations (this result for Cair) 
can be derived in a more controlled way using the Exact Renormalization Group |46|). 
Concerning chaos properties, it was shown in Ref. [44| that Gi2(0) grows linearly with 
a for small £ before it saturates to a constant for large £, Gi2(0) ~ > 0, which yields 
Si2{q) oc a/q^. These results close to Tg can be summarized as 

{logl/g 1 
^— ^ , 9 » ^5 
a ^^^^ 
^ ' ^ « ^s' 

while in real space, the behavior of the interlayer correlation function Ci2{r) for T <Tg 
is thus 

C,.(r)~|"'°«''''' • ''^^^^ (24) 
[(Tlog(r) , r Ls . 

At T = 0, it was recently shown [34], using FRG to one loop including the term 
Gii(m) that the intralayer correlation function also behaves like C(r) oc log^(r), in rather 
good agreement with numerics [HI [32]. One thus also expects that Ci2(r) oc log^(r) 
for r < L5 [m [321 [34]. For r > L5, a behavior of C^ir) ~ alogr as in Eq. ((231) 
was discussed in Ref. [43]. To determine analytically whether a > at T = requires 
a detailed and difficult analysis of the coupled FRG equations for Rij{u),Gij{u) ([22l) . 
which goes beyond the previous studies done in that direction in Ref. [3H[43l[45]. Here, 
we will answer this question using numerical simulations. 

The purpose of our study is actually to answer the two main questions : (i) what 
are the residual correlations beyond Ls and in particular is a also finite at T = ? (ii) 
what is the scaling form of this correlation function Ci2(r), i.e. the analogous of Eq. 
( flTl) from which one can extract the overlap length Ls and check the value of chaos 
exponent a = 1 as expected from droplet scaling ? 

Here will use the minimum-cost-fiow algorithm described above to compute the 
two ground states h} and hf with free boundary conditions. Instead of the correlation 
function Ci2{r) we compute numerically the Fourier transform S^iq) of the "overlap" 
between the configurations (fTGll . Our simulations have been performed on a square 
lattice of linear size L = 256 and we have chosen q = (g, 0) with q = 27in/L with 
n = 0,l,2,---,L — 1. The disorder average has been performed over 10^ independent 
samples. In Fig. [H] (left) we show our numerical data for Sui^q). Its behavior close 
to Tg ([23l) suggests strongly to plot q^Si2{q) as a function of q. Given that we are 
working on a discrete lattice of finite size, it is more convenient to work with the 
variable y = sin (g/2) = [(1 - cos (g))/2]i/2 instead of q (of course for small q it makes 
no difference). In Fig. [IT] (left) we actually show a plot of S u (q) [sin {q / 2)]"^ /A{6) as a 
function of sin(g/2) on a log-linear plot for different values of 5 = 0,0.1,0.2,0.3. On 
this plot the amplitude A{6) is chosen such that the curves for different values of 6 do 
coincide for sin(g/2) ~ 1, with ^4(0) = 1. For 5 = 0, where Suiq) = Su{q), this plot 
is almost a straight line, which suggests indeed that for small q where sin (g/2) ~ g/2. 
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Figure 11. Left : Plot of S 12 {q) * y'^ /A{S) as a function oiy = sin (q/2) for different 
values of (5 = 0, 0.1, 0.2, 0.3 for a linear system size L = 256 on a log-linear scale. The 
deviation from the straight line for (5 > is a clear indication of disorder chaos. This 
behavior is consistent with the one in Eq. l(23|l and in particular with a finite value of 
(T at T = 0. Right : Same data as the one shown in the right panel where we plot 
v'^ Si2{q) / A{5) + alogS as a function of y/5. The relatively good collapse is consitent 
with the scaling form proposed in Eq. l(25|) . implying in particular a = 1. 



Siiio) ~ log (l/?')/?'^- this produces the log^(r) behavior of the correlation function 
Cii(r) in real space. This result is consistent with previous numerical studies of the 
ground state 

For (5 > the behavior is however quite different : indeed for small g, the curves 
deviate from the straight line indicating a saturation to a finite value. We have checked 
that this is not a finite size effect: in particular, for a given value of 5, we have checked 
that the value of q where the bending occurs and the saturating value do not depend on 
L. This bending is thus a clear signature of disorder chaos in this model. In addition, 
the saturation to a finite value is consistent value of a at T = 0. If one translates these 
results in real space, one obtains a behavior of the correlation function Ci2(r) oc alogr 
as in Eq. f l24ll . Of course A{5) — > 1 when but we were not able to characterize 

precisely the dependence of A{5). We emphasize that the observation of disorder chaos 
here needs a precise computation of q^Si2{q) and thus needs accurate statistics. 

The scaling form in Eq. (fTSll would suggest to plot q^Si2{q) as a function of 
L'^'^(p{qLs). Here one has C = together with a = 1 expected from scaling argument 
and on the other hand one expects, see Fig. [TT] (left) and also Eq. (l23l) . that 
(p{x) ~ — alog(x) for large x (with a = 0.016(1) estimated from q'^Sii{q)). Therefore, 
to guarantee that Sui^q) has a good limit when 5 — 0, we propose the scaling form 

5i2(g) ~ g-' (^) - a log 5) , (25) 

§ In addition, this gives an alternative way to estimate the amplitude of the log^(r) term, here 02 — 0.40 
which is sHghtly different from the previous estimate obtained from the fit of Cii(r) in real space and 
yielded a2 = 0.57 (which actually recovered in the present simulations by applying the same fitting 
procedure). 
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where the scaling function ip{z) behaves like 




—a log z , z ^ oo 



(26) 



such that the additional constant —a log 5 with a = 0.016(1) (independently of 6) is 
needed to yield a well defined limit 5 — > 0. We have checked this scaling form (1251 ) 
for different values of 6 = 0.1,0.2,0.3. Again, to take into account finite size effects, 
we use the variable y = sin(g/2) instead of q. In Fig. [TT] (right) we show a plot of 
y"^ Si2{q) / A{6) + a log 5 as a function of y/6. The relative good collapse of these curves 
corresponding to different values of 6 is in rather good agreement with the scaling in 
Eq. (I25ll . This indicates in particular that a ~ 1, although a study of Suiq) for smaller 
values 6 would certainly be necessary to obtain a more precise estimate of this exponent. 

5. Discussion 

We found that the left passage probability of zero-temperature domain walls in the 
disordered SOS model in different geometries agrees well (within the numerical error 
bars) with Schramm's formula for k = 4. Since the fractal dimension of the domain 
walls is ds = 1.25 ± 0.01 this is inconsistent with ds = I + k/8 = 1.5 that should hold 
if this ensemble of random curves would be described by SLE. One condition for SLE 
to hold is that the measure of the ensemble of random curves fulfills a domain Markov 
property. To check it numerically is a computationally extremely challenging endeavour 
(c.f. |23j) and we did not attempt it here. The second condition for SLE is that the 
measure is invariant under a conformal mapping of the domain within which the random 
curves are defined. In our study we considered several different domains shapes and we 
found that in all cases the left passage probability is described by Schramm's formula 
adapted via a conformal map to the specific geometry under consideration. Although not 
being a sufficient criterion for conformal invariance to hold in general this observation 
is at least surprising regarding the fact that k and ds do not obey the SLE prediction. 
To shed light on this issue further tests of conformal invariance are worthwhile. 

Interestingly the contour-loops in this system (i.e. the lines connecting sites of equal 
height in the ground state) have a fractal dimension close to df = 1.5 |39j, the same 
dimension as the contour loops in the SOS model without disorder at finite temperatures 
above the critical temperature [22l [40]. The latter are indeed described by SLE with 
K = 4 [22], [40]. It remains to be checked, whether the contour loops in the disordered 
SOS model are also also described by SLE with k = 4. 

The boundaries of droplets, i.e. connected clusters of a given lateral size that have a 
minimal excitation energy, have the same fractal dimension as domain walls {dg = 1-25), 
but their energy AEi saturates at a small, finite value for increasing lateral size /. In 
[20] it was shown that the mean droplet energy AE^ decreases with system size L like 
AE^^ ~ InL. This result is actually consistent with ours, since in [20] the droplet 
size was only restricted by an upper bound (the system size L) but not by a lower 
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bound as in our study. So the result for AE^ is the optimal excitation among one of 
typically \ogL excitations of fixed (maximum and minimum) size /. //the distribution 
of the energies of these excitations on different length scales is independent of the length 
scale, then AEl is just the minimum of InL independent, identically distributed random 
numbers, thus proportional to 1/ InL. What we show in our study is essentially that this 
assumption is indeed fulfilled, as shown in Fig. [TOl This is a non-trivial result and does 
not immediately follow from a vanishing stiffness exponent 6 = 0. With a view towards 
their entropic contribution one would like to know how the number of independent (i.e. 
spatially disjoint) droplets scales with their size. It should be possible to study these 
challenging questoin with the methods that we presened. 

Finally, we have shown that there is disorder chaos at T = in this model. Our 
numerical data show a behavior of interlayer correlations S^iq) compatible with Eq. 
( l23ll with (T > 0, thus rather similar to the one obtained close to Tg [44]. Here, the two 
ground states of the system, embedded in two slightly different disorder realisations, 
thus display logarithmic residual correlations. In addition, our data are consistent with 
a chaos exponent a = 1, in agreement with droplet scaling argument. Following the 
studies of Ref. [34l[43l[45] it would certainly be interesting to describe analytically these 
behaviors by analysing in detail the model in Eq. (l22ll at T = 0. Of course, especially in 
view of recent analytical progress [45], it would also be very interesting to extend these 
studies of disorder chaos at finite temperature, which we expect to play a role here in 
this marginal glass phase. 
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